AD-A251  087 


REPOR 


Form  Approved 
0MB  NO.  0704-0] 88 


’•joiit  -eocrtinc  5L/3fn  <0'  t.ii*  colicc  nq  th«  time  lOf  reviewmq  instruaiont.  teerchinq  onlinq  am  iourcn. 

'iinerirq  jne  .^einieminq  me  dei>  n«  d  comments  reqerdinq  this  burden  estimate  or  any  other  asoect  a<  inis 

collection  of  ini’rmation.  inciudinq  su  ts.  Directorate  tor  information  Operations  and  heoorts.  i  i  15  lelferson 

Dans  Hiqhwa*.  Suite  IJoi.  Arlinqton.  V  •»  raaur->]ur.  and  to  the  OHiceoi  Manaqement  and  tudqet.  Papcrworii  Reduction  Project  (0704-0  t8S).Washinqton.  DC  7050]. 


1.  AGENCY  USE  ONLY  (Leive  bUnk) 


4.  TITLE  AND  SUBTITLE 


2.  REPORT  DATE 

May  1992 


3.  REPORT  TYPE  AND  OATES  COVERED 

Technical  Rep.  (Interim)  7/1/91-6/30/92 


5.  PUNOING  NUMBERS 


Dynsunics  of  Ion  Solvation  in  a  Stockmayer  Fluid 


Grant  N00014-89-J-1169 
RT  4131036 


6.  AUTHOR(S} 

L.  Perera  and  M.  Berkowitz 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORE 

University  of  North  Carolina 
Department  of  Chemistry 
CB#  3290 
Chapel  Hill,  nc  27599 


DTIC 


JUN3  1992 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

Technical  Report  #7 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 

Office  of  Naval  Research 
800  N . fQuincy  Street 
Arlington,  VA  22217-5000 
(Code  1113PS,  P.P.  Schmidt) 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 


Published  in  Journal  of  Chemical  Physics 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


This  document  has  been  approved  for  public  release  and 
sale;  its  distribution  is  unlimited. 


14.  SUB/ECT  TERMS 


IS.  NUMBER  OF  PAGES 


Stockmayer  fluid,  ion  solvation,  solvent  relaxation  time 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  ‘19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

Unclassified  Unclassified  Unclassified 


Unclassified 


NSN  7S40-01-280-$S00 


Standard  Form  298  <R*v  2-89) 
bmcnMd  by  aNV  5M.  tn-'l 


OFFICE  OF  NAVAL  RESEARCH 


Grant  N00014-89-J-1169 
R&T  Code  4131036 
Technical  Report  #  7 


Ac«««3io«  For 
NT  IS  Gilifcj: 

tiii 

Ju:.  t  if  laa*  1  on. 


□ 

□ 


!  iv. ... 

1  Dlstr 

;  Avsilabtlity  Codes  i 

Dlst 

AtbII  and/or 
Special 

”  Dynamics  of  Ion  Solvation  in  a  Stockmayer  Fluid.” 

by 

L.  Perera  and  Max  L.  Berkowitz 


Published  in 

Journal  of  Chemical  Phsrsics  3092  (1992) 


The  University  of  North  Carolina  at  Chapel  Hill 
Department  of  Chemistry 
Chapel  mu,  NC  27599-3290 


MAY  1992 


Reproduction  in  whole  or  in  part  is  permitted  for  any  purpose  of  the  United  States 
Government 

This  docmnent  has  been  approved  for  pubUc  release  and  sale;  its  distribution  is  unlimited. 


92-14403 

liiilllill 


92  6  01  0  72 


1 


October  1991  -t-  _ 

Jjt  ore^s 

J.  Chem.  Phys.  (Submitted) ' 


DYNAMICS  OF  ION  SOLVATION  IN  A  STOCKMAYER  FLUID 
Lalith  Perera  and  Max  L  Berkowitz 

Department  of  Chemistry,  University  of  North  Carolina,  Chapel  Hill,  North  Carolina,  27599 


Abstract:  Molecular  dynamics  computer  simulations  were  performed  to  study  the  dynamics 
of  the  ionic  solvation  in  a  Stockmayer  fluid.  The  simulations  show  that  the  solvent  relaxation 
proceeds  in  two  time  regimes.  Most  of  the  relaxation  occurs  in  a  short  time  period  during 
which  the  relaxation  process  can  be  described  a  gaussian  function.  The  long  time  regime 
can  be  described  by  an  exponential  relaxation.  The  decay  exponent  of  the  relaxation  fimction 
in  this  regime  is  the  same  as  the  exponent  describing  the  decay  of  the  single  dipole 
correlation  function.  In  addition,  the  contribution  of  the  rotational  and  translational  modes 
of  the  solvent  to  the  energy  relaxation  was  investigated.  It  was  found  that  when  the 
rotational  mode  is  the  dominant  mode  of  the  solvent  motion  the  relaxation  occurs  from  the 
outside-in,  in  accordance  with  the  Onsager  "snowball”  picture.  When  the  influence  of  the 
translational  mode  is  increased  the  Onsager  picture  breaks  down. 
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Introduction 

The  dynamic  response  of  the  solvent  to  the  sudden  change  of  the  charge  distribution 
in  a  solute  molecule  has  been  the  subject  of  many  investigations,  theoretical  and  experimental 
[1].  While  the  experimental  data  are  obtained  from  the  applications  of  &st  laser  pulse 
techniques  to  rather  complicated  systems  [2,3],  the  majority  of  analytical  theories  use  rather 
simplified  models  to  describe  the  solutions  [4-10].  At  the  same  time  the  computer  simulations 
are  performed  on  models  that  use  a  simple  description  of  the  solute  and  a  rather  realistic 
description  of  the  solvent  [11-16].  Thus,  re<»nt  simulations  were  performed  to  study  the 
response  in  solvents  such  as  water[ll,12],  methanol  [13]  and  acetonitrile  [14].  To  compare 
the  results  of  the  simulations  with  the  analytical  type  theories,  it  is  desirable  to  perform  a 
simulation,  where  both  solute  and  solvent  are  described  by  simple  potentials.  Therefore,  in 
this  paper  we  present  results  frnm  such  a  simulation  and  compare  them  with  the  results  from 
the  existing  analytical  theories,  such  as  the  dynamical  mean-spherical  approximation  (MSA) 
[6-8]  and  theories  based  on  the  Smoluchowski-Vlasov  Equation  (SVE)  [8-10].  We  also  address 
the  issue  of  the  validity  of  Onsager’s  hypothecs  [16],  and  the  role  of  the  translational  and 
rotational  motions  in  the  relaxation.  Finally,  we  express  the  hope,  that  simulations  performed 
on  simple  systems  can  serve  as  a  reference  point  for  the  development  of  more  sophisticated 
molecular  theories,  such  as  the  one  given  in  ref  [17]. 

Model  and  Methodology 

Our  system  consists  of  one  solute  particle  dissolved  in  a  Stockmayer  solvent  [18].  The 
choice  of  this  particular  model  is  primarily  due  to  the  £Eu:t  that  analytical  theories  based  on 
the  dynamical  extension  of  the  MSA  model  [6-8]  and  on  the  SVE  [8-10]  are  available  for  a 
model  very  similar  to  ours  (the  theories  are  done  for  hard  spheres  with  embedded  dipoles). 
In  addition  the  detailed  knowledge  of  the  dielectric  properties  of  thia  solvent  is  available  from 
molecular  d3mamic8  simulations  [19-21]. 

The  interaction  energy  for  a  pair  of  Stockmayer  particles  is  given  by  the  following 
expression: 


Uffir)  -  4«((a/r)“-(a/r)®]  + 


(1) 
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where  r  is  the  distance  between  particles  i  and  j,  Mt  is  the  dipole  moment  of  particle  i,  is 
the  dipole  moment  of  particle  j,  e  and  a  are  the  Lennard-Jones  parameters.  The  dipole 
tensor  Ty  has  its  usual  form: 

T,  - 


where  I  is  the  unit  tensor. 

The  solute-solvent  interaction  is  represented  by  the  term: 

U/r,nj)  -  46C(<y/r)“-(a/r)®]-#tjt®  (3) 


where  E  is  the  electric  field  due  to  the  solute  charge,  and  is  given  by  an  equation  E  »  qr^/ry^. 
For  the  neutral  solute,  this  quantity  becomes  zero  and  the  last  term  in  the  equation  (3) 
vanishes.  The  equations  of  motion  for  our  q^stem  are  obtained  firom  the  Tiagrangian 


»>/  j 


(4) 

where  the  first  two  terms  represent  translational  and  rotational  kinetic  energies  of  the  solvent 
molecules,  the  third  term  represents  the  dipole-dipole  interaction  energy  of  the  solvent,  the 
fourth  term  is  representing  the  Lennard-Jones  interaction  energy  of  the  system.  To  preserve 
the  magnitude  of  the  dipole  moment  of  the  solvent  we  use  the  method  of  constraints  [22]. 
Therefore  the  fifth  term  of  the  Lagrangian  contains  a  constraint  variable  with  the 
corresponding  Lagrange  multiplier.  The  last  term  in  the  Lagrangian  is  due  to  the  solute- 
solvent  electrostatic  interaction.  The  solvent  is  characterized  ( in  reduced  units)  by  a  dipole 
moment  >3.0,  a  moment  of  inertia  V  >0.026,  a  density  p*  >0.822  and  a  temperature  T* 

>  1.15.  The  Lennard-Jones  parameters  for  solute-solvent  interaction  were  taken  to  be  the 
aama  as  for  solvent-solvent  interaction.  The  simulations  were  carried  out  in  a  cubic  box, 
which  contained  511  solvent  molecules  and  one  solute  molecule.  For  simplification  of  the 
treatment,  the  solute  particle  was  considered  to  be  immobile.  Equations  of  motion  for  solvent 
molecules  were  integrated  using  the  leapfing  algorithm  and  the  rotational  part  of  the  motion 
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was  handled  by  the  methods  of  constraints.  A  good  description  on  the  implementation  of  this 
method  is  found  in  Ref.  [23].  The  step  size  in  reduced  units  of  time  was  chosen  to  be  equal 
to  0.0025.  To  facilitate  the  comparison  of  our  results  with  the  results  from  other  simulations 
note  that  this  time  step  is  equal  to  5.4  fs  if  the  argon  parameters  (e  » 119.8K  and  a  *  3.405A) 
are  chosen  for  solvent  particles.  Periodic  boundary  conditions  along  with  the  spherical  cut-off 
convention  were  implemented  in  the  calculations  and  the  reaction  field  technique  was  used 
to  account  for  the  long  range  dipole-dipole  interactions.  In  the  reaction  field  calculations  we 
used  the  value  of  the  bulk  dielectric  constant  of  the  Stockmayer  fluid,  which  is  66.1  for  the 
values  of  the  reduced  parameters  given  above  [20]. 

In  the  present  work,  we  carried  out  three  equilibrium  simulations  and  three  sets  of 
non-equilibrium  simulations.  The  first  equilibrirun  simulation  (ESI)  was  performed  for  the 
^tem  with  the  uncharged  solute,  while  the  solute  was  charged  (q>le)  in  the  second 
simulation  (ESII).  £a  the  first  two  simulations  the  reduced  moment  of  inertia  T  of  the  solvent 
molecule  was  assigned  the  value  of  0.025,  in  the  third  equilibrium  simulation  the  value  of  I* 
was  increased.  By  increasing  the  moment  of  inertia  one  achieves  slow  rotational  motion  of 
the  solvent  and  thus  amplifies  the  effects  of  translational  motion  on  the  overall  dynamics. 
Since  we  do  not  intend  to  mimic  the  behavior  of  any  real  solvent  S3rstem,  we  can  make 
arbitrary  changes  in  some  parameters  to  gain  more  physical  insight  into  the  processes  w? 
discuss  here.  The  third  equilibrium  simulation  (ESIQ)  was  performed  using  mostly  the  samp 
set  of  parameters  as  used  in  ESI,  only  I*  was  increased  100  times.  In  all  the  equilibrium 
simulations  the  trajectories  consisted  of  20  ps  equilibration  period  followed  by  100  ps 
production  run.  The  average  reduced  temperature  of  each  trajectory  was  maintained  at 
T* « 1.15  by  occasional  rescaling  of  the  velocities. 

To  mimic  the  solvation  phenomena  we  use  the  non-equilibrium  simulations 
technique.  In  these  simulations  the  solute  chai^  distribution  is  instantaneously  changed  and 
the  response  of  the  solute  is  measured.  To  quantify  the  solvation  dynamics  one  calculates  in 
the  simulation  the  normalized  response  fimction  or  relaxation  fimction 


_  SEit)— SE(po) 
"  SEm-SE(co) 


(5) 
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In  equation  (5)  £E(t)  is  the  fluctuation  in  the  solvent-solute  interaction  energy  at  a  given  time 
t,  averaged  over  all  the  trajectories  in  a  particular  set  of  simulations.  The  equilibrium 
trajectories  of  the  system  with  the  neutral  solute  are  used  to  generate  initial  configurations 
for  the  non-equilibrium  trajectories.  Therefore,  we  have  selected  50  configurations  (separated 
by  2  ps  each)  firom  the  equilibrium  simulation  SEI  and  at  t«0  we  placed  a  charge  of  +  le  on 
the  solute.  To  compare  the  results  of  our  simulations  with  some  analytical  theories,  such  as 
MSA,  in  which  the  relaxational  mechanism  is  due  to  the  rotational  motion  only,  we  fi-oze  the 
translational  motion  of  the  molecules  in  the  first  set  of  50  trajectories.  This  first  set  of  50 
non-equilibritun  tngectories  we  call  NESI  (non-equilibrium  set  L)  To  fiieeze  the  translations 
we  imposed  very  heavy  masses  on  the  solvent  molecules  (10^°  times  the  actual  mass).  Also, 
the  initial  linear  velodties  were  set  equal  to  zero  to  ensure  that  particles  do  not  move 
initially.  After  the  instantaneous  charge  jump,  the  tnyectories  fi:t>m  the  molecular  dynamics 
were  monitored  for  5.4  ps.  (The  equilibration  has  been  achieved  within  2.0  ps). 

It  has  been  emphasized  that  translational  motion  of  the  solvent  molecules  plays  a 
substantial  role  in  the  overall  mechanism  of  the  relaxation.  This  was  confirmed  recently  by 
Bagchi  and  Chandra  [10],  who  extended  the  SVE  treatment  to  indude  the  effects  of  the 
translational  motioiL  Thqy  demonstrated  that  if  translational  diffhsion  is  considerably  large, 
the  Onsager  picture  of  the  solvent  relaxation  occurring  in  a  'snowball  ”  fiuhion  is  no  longer 
acceptable  [24].  We  tried  to  address  thi«  issue  in  two  other  sets  of  non-equilibrimn 
simulations.  In  one  set  of  simulations  (NESff),  consisting  again  of  50  trajectories,  we  started 
the  trajectories  firom  the  same  initial  positions  as  in  set  NESL  The  initial  velocities  of  the 
solvent  molecules,  which  were  set  equal  to  zero  in  the  NESI,  were  now  assigned  fi-om  the 
equilibrium  trtuectories.  The  lengths  of  the  trajectories  were  restricted  to  2.7  ps  since  the 
equilibration  was  attained  within  2  pe. 

The  third  set  of  non-equilibrium  trajectories  (NESm)  was  started  firom  50  points  of 
ESm  trajectory,  where  the  moment  of  inertia  of  the  solvent  was  100  times  larger  than  that 
used  in  ESL  Trajectories  were  continued  for  22  ps  after  the  initial  charge  jump.  No  attempt 
has  been  made  to  control  the  temperature  of  the  ^stem  in  any  of  the  above  mentioned  non¬ 
equilibrium  tngectories. 
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Results  and  Discussion 

Using  the  data  &om  equilibrium  simulations,  we  first  consider  the  structural  properties 
of  solvent  around  the  solute  molecule.  The  solute-solvent  radial  distribution  fimctions  (rdf  s) 
obtained  firom  ESI  and  ESU  are  shown  in  Figure  L  The  curves  on  this  figure  demonstrate 
the  existence  of  the  well  defined  solvation  shell  structtnres  around  the  solute,  regardless  of 
the  charge  on  the  solute.  In  Table  1  we  summarize  the  information  extracted  fi^om  the  rdf  s 
along  with  the  other  relevant  data.  As  the  rdf  s  in  Figure  1  indicate  we  distinguish  three 
solvation  shells  around  the  solute.  The  rest  of  the  solvent  is  ”bulk-like”.  The  boundaries  of 
the  shells  are  given  by  the  locations  of  the  minima  of  the  rdfk  The  properties  shown  in  Table 
1  are  averages  over  the  shells  defined  by  such  boundaries.  The  orientatioxial  distribution  of 
the  angle  between  the  solvent  dipole  and  a  vector  connecting  the  center  of  a  solute  and  the 
center  of  the  solvent  is  shown  in  figure  2. 

From  Figures  1  and  2  and  Table  1  we  can  see  that  the  incorporation  of  the  charge  on 
the  solute  site  results  in  the  increase  of  the  sharpness  of  the  first  peak  of  the  rdf,  while  the 
peak  position  shifts  inward.  The  size  of  the  first  shell  decreases  and  as  a  result  the  number 
of  nearest  neighbors  decreases.  The  same  is  true  for  the  second  sheU.  At  the  same  time, 
note,  that  the  denaty  of  solvent  in  the  first  shell  around  the  ion  is  larger  than  the  density 
of  the  solvent  in  the  first  shell  around  the  neutral  solute.  The  orientational  distribution  is 
homogeneous  in  the  case  of  solvated  neutral  solute  but  displays  a  strong  orientational 
preference  in  the  case  when  the  solute  is  charged.  As  we  can  see  in  this  case  the 
orientational  preference  holds  for  all  3  shells  around  the  solute.  The  rdf  s  and  orientational 
distribution  fimctions  given  by  Figures  1  and  2  represent  the  iTiitial  and  final  states  of  our 
non-equilibrium  simulations. 

As  we  have  already  mentioned,  the  response  of  the  solvent  to  the  instantaneous 
change  of  the  solute  charge  distribution  is  measured  in  a  molecular  dynamics  computer 
experiment  by  the  normalized  response  fimction  (see  equation  5). 

In  Figure  3  we  display  the  results  for  the  correlation  fimction  S(t)  obtained  fi-om  the 
NESI  simulation  along  with  the  results  obtained  fi’om  dynamical  MSA  theory  [7]  and  SVE 
type  theory  [10].  The  functional  forms  of  the  relaxation  fimction  S(t)  obtained  in  the 
fiamework  of  MSA  and  SVE  and  used  here  are  presented  in  the  Appendix.  Since  MSA  theory 
does  not  take  into  account  the  translational  mode  of  the  relaxation,  the  comparison  of  this 
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theory  with  the  data  from  NESI  simulation  is  most  appropriate.  The  SVER  curve  was 
obtained  by  solving  SVE  where  only  the  rotational  mode  of  the  relaxation  was  included. 
Finally,  to  complete  the  picture,  we  have  shown  in  the  figure  the  exponential  relaxation  of 
the  solvent  predicted  by  the  continuum  theory. 

As  we  observe  in  Figure  3  the  relaxation  curve  corresponding  to  NESI  displajrs  a  fast 
initial  decay  (~0.1ps)  followed  a  relatively  slow  decay  ('>2.0ps)  modulated  by  periodic 
oscillations.  We  also  observe  that  after  the  initial  response  of  the  solvent  the  energy  of  the 
^tem  is  even  lower  than  its  equilibrium  value  (Le.  S(t)  becomes  negative).  This  may  be 
attributed  to  a  quick  initial  reorientation  of  dipoles  towards  the  ion,  which  results  in  the 
lowering  of  the  ion-dipole  energy.  But  the  system  does  not  remain  in  this  configuration  with 
the  lowest  ion-dipole  energy  due  to  the  presence  of  dipole-dipole  interactions.  Therefore,  the 
observed  oscillations  in  the  correlation  function  are  the  results  of  the  timely  adjustment  of 
the  dipoles  to  the  field  of  the  fireshly  crMted  ion  and  the  field  emanating  from  the 
neighboring  dipoles. 

Let  us  now  compare  the  relaxation  behavior  predicted  by  the  continuum,  MSA  and 
SVE  theories  with  the  molecular  dynamics  results.  Perhaps  we  shall  mention  here,  that  our 
solvent  is  quite  accurately  described  by  the  continuiun  Debye  model  [19].  The  simplest  model 
for  the  relaxation,  the  model  in  which  the  solvent  is  assumed  to  be  represented  by  a 
continuum,  predicts  an  exponential  decaying  relaxation  with  a  time  constant  corresponding 
to  the  longitudinal  relaxation  time  r,  >  Substitution  into  this  formula  of  the 

values  2.15p8,  66.1  and  LO  of  the  Debye  relaxation  time,  staticfcg)  and  opticaKe  J  dielectric 
constants  respectively  [20],  predicts  in  our  case  a  relatively  short  relaxation  time  (0.033ps). 
As  Rips  et  aL  pointed  out,  this  time  provides  the  lower  limit  to  the  observed  relaxation  [7]. 
However,  it  is  necessary  to  mention  here  that  for  most  of  the  experimental  systems  this  time 
provides  a  reasonably  good  estimate  to  the  relaxation  time  [3]. 

A  more  sophisticated  theory  that  takes  the  solvent  structure  into  accoimt  is  the 
dynamical  MSA  theory  [7].  It  predicts  a  multiezponential  relaxation  behavior  with  relaxation 
times  ranging  between  r,  and  r^.  Although  the  dynamical  MSA  theory  considers  the 
molecular  character  of  the  solvent,  it  does  it  in  a  very  average  way  and  as  a  result  no 
osdllatoiy  behavior  of  the  relaxation  is  predicted  by  it.  Also,  as  Figure  3  shows,  the  initial 
behavior  of  S(t)  predicted  by  MSA  is  incorrect  From  the  MD  simulations  we  condude  that 
the  initial  decay  may  be  better  characterized  by  a  gaussian  function  than  by  an  exponential 
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one.  Such  a  gaumrian  decay  of  the  initial  response  has  been  observed  in  other  simulation 
works  [15].  In  addition  we  also  observe  from  Figure  3,  that  the  time  to  approach  the 
asymptotic  limit  (S(t)»0)  is  much  longer  in  the  framework  of  MSA  than  what  MD  shows. 

The  relaxation  curve  obtained  firom  the  theory  based  on  the  use  of  SVE  [10]  is  close 
in  our  case  to  the  predictions  of  the  simple  continuum  modeL  As  we  can  see,  the  curve  from 
the  molecular  dynamics  displays  a  different  behavior  at  the  initial  time  but  its  aqnnptotic 
behavior  is  the  w«mA  as  predicted  by  models  based  on  SVE  and  continuum  description.  In 
general  we  can  say  the  molecular  djniamics  displays  a  behavior  somewhere  between  that 

predicted  by  the  continuum  theory,  MSA  and  SVE.  If  one  describes  the  relaxation  by  a  single 
relaxation  time  r’  defined  by: 


r 


(6) 


one  obtains  firom  the  molecular  (fynamics  sirmilation  that  r*  is  around  0.2p8.  The  MSA 
treatment  predicts  a  value  of  0.17  ps,  which  is  in  reasonable  agreement  with  the  relaxation 
time  obtained  from  the  molecular  dynamics.  But  this  agreement  should  be  considered  as 
fortuitous,  no  special  physical  meaning  should  be  ascribed  to  it. 

Let  us  now  study  how  the  inclusion  of  solvent  translational  motion  affects  its 
relaxation  dynamics.  In  Figure  4.a  we  compare  the  results  obtained  from  a  simulation  in 
which  the  translational  motion  is  included  in  the  dynamics  (NESIl)  with  that  of  the  previous 
sim  'lation  (NESI).  Except  for  certain  details  discussed  below,  the  two  curves  display  similar 
behavior.  This  leads  us  to  a  conclusion  that  for  this  particular  system,  the  rotational  motion 
of  the  solvent  essentially  governs  the  overall  relaxation  mechanism.  At  least  for  this  system, 
it  is  not  surprising  to  observe  such  behavior,  since  it  has  a  relatively  large  rotational  diffrision 
constant  (D,  «  L167/ps)  compared  to  its  translational  diffusion  constant  (D^  «  0.42A^/ps). 

What  is  the  origin  of  the  different  behavior  of  curves  NESI  and  NESn  in  Figure  4.a, 
especially  in  the  region  around  the  first  minimum?  We  have  observed  finm  the  dynamics  that 
the  reconstruction  of  the  solvent  shells  around  the  ion  (the  electrostriction  effect)  takes  place 
within  the  first  0.2ps  at  which  the  ioitial  fast  relaxation  occurs.  The  reconstruction  of  the 
shell  structure  brings  the  molecules  in  the  first  two  shells  much  closer  and  thus,  causes  the 
solvent  molecules  to  feel  the  repulsions  fiiom  their  neighbors  somewhat  stronger  than  in  the 
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case  when  the  translation  is  frozen.  The  rea>nstruction  results  in  an  earlier  activation  and 
stronger  dipolar  interactions  between  solvent  molecules.  This  is  especially  visible  from  the 
curves  in  the  region  around  O.lSps.  However,  the  reconstruction  does  not  seem  to  affect  the 
long  time  relaxation  behavior. 

Figure  4.a  also  includes  a  curve  labeled  SVERT.  This  curve  was  obtained  from  the 
solution  of  SVE  [10]  where  the  effect  of  rotational  and  translational  diffusion  were  included. 
Again  only  long  time  behavior  of  the  SVERT  curve  is  correct,  as  expected.  To  study  the 
relative  contributions  of  translational  or  rotational  diffusional  modes  into  the  total  relaxation 
process  a  parameter  p’  defined  by  an  equation 

p'  -  DJ(2D^  (7) 

was  introduced  [10].  Our  Stockmayer  solvent  in  which  T »0.025  yields  a  value  of  p’»0.012. 
As  was  pointed  out  Bagchi  and  Chandra  this  value  of  p’  should  not  change  the  relative 
importance  of  translational  and  rotational  modes  [10].  To  study  how  the  translational  motion 
of  the  sohrexit  influences  the  relaxation  dynamics  at  larger  values  of  p’,  we  suppressed  the  last 
rotations  of  solvent  molecules  by  increasing  their  moment  of  inertia  by  a  factor  of  100  and 
performed  an  equilibrium  (ESllll  and  a  series  of  non-equilibrium  simulations  (NESm)  for  this 
new  system.  The  equilibrium  simulations  show  that  the  translational  diffusion  coefficient 
remains  near  that  of  the  original  system  (0.38AVps),  but  the  rotational  diffusion  coefficient 
decreases  by  a  factor  of  38  ( to  the  value  0.031/ps).  For  this  new  Stockmayer  solvent  p’  >  0.42, 
and  therefore  we  expect  to  see  a  more  pronounced  effect  of  translation  on  the  relaxation.  If 
we  would  fiieeze  out  completely  the  translational  motion  of  the  solvent,  the  dynamics  in 
NESm  would  be  scaled  up  by  a  factor  of  10.  Consequently,  one  would  expect  to  see  a 
complete  relaxation  of  the  solvation  around  the  ion  to  take  place  after  around  17-18  ps. 
However,  as  Figure  4.b  indicates,  a  complete  relaxation  is  achieved  after  around  7  ps.  This 
points  to  an  appreciable  contribution  of  translational  motion  to  the  total  relaxation  process. 
Another  effect  that  tranalational  motion  has  on  the  relaxation  is  related  to  the  Onsager 
"snowball*  cozyecture  and  ia  discussed  below. 

As  we  hove  seen  firam  our  simulation  and  other  simulations,  the  relaxation  of  the 
system  can  be  described  by  a  fast  initial  decay  and  a  subsequent  slow  relaxation.  The 
existence  of  these  two  time  scales  are  also  confirmed  by  experiments  [25].  Therefore  the 
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main  questions  we  should  ask  are:  what  are  these  two  (or  more)  time  scales  observed  in 
solvation  dynamics  and  how  do  they  correlate  with  the  dynamics  of  solvent  molecules?  In  the 
previous  simulations,  the  short  time  decay  is  assigned  to  the  inertial  motion  of  the  solvent 
molecules  [11,13,14],  meaning  that  the  motion  of  the  solvent  can  not  be  described  by  the 
diffusion  equation.  Indeed  we  observe  that  the  initial  fast  decay  is  essentially  due  to  the 
immediate  reaction  of  the  solvent  molecules,  especially  of  the  first  solvation  layer.  The 
nature  of  the  long  time  decay  is  clarified  by  Figure  5,  where  the  correlation  function  S(t)  is 
compared  to  the  single  dipole  autocorrelation  function  4(t),  defined  as: 

#(f)  -  <M(t)*M(0)>  (8) 

</i*> 


The  autocorrelation  function  4(t)  was  calculated  from  a  separate  simulation  performed  on  512 
Stockmayer  particles.  As  we  can  see  from  Figure  5  the  long  time  decay  of  S(t)  is 
characterized  by  the  same  relaxation  time  as  the  decay  of  4(t).  That  does  not  mean  that  the 
long'time  dynamics  of  a  dipole  is  determined  by  solving  the  dynamical  equations  for  a  single 
particle,  since  to  find  a  sint^e  particle  correlation  function  in  the  many>body  problem  one  has 
to  solve  the  naany-bodty  problem.  The  agreement  in  long-time  behavior  of  S(t)  and  4(t) 
provides  a  reasonahle  conformation  to  the  idea  that  the  long  time  decay  of  the  relaxation  is 
due  to  the  a4justment  of  the  dipole  orientation  to  the  resulting  field  from  the  rest  of  the 
S3rstem. 

At  this  stage  we  want  to  return  to  the  discussion  of  the  Onsager  "snowball"  effect. 
Onsager  once  commented  that  the  relaxation  of  solvent  around  a  newly  created  charge 
distribution  will  proceed  from  outside  towards  the  solute  -  like  a  "snowball"  [16].  We  can  check 
Onsager’s  conjecture,  since  molecular  dynamics  simulations  allow  us  to  examine  the 
contributions  to  the  relaxation  emanating  from  the  different  solvent  shells  aroimd  the  ion. 
Thus  in  Figure  6a  we  show  the  plots  of  S(t)  vs  time  for  each  solvent  shell  for  NESL  Initially, 
all  the  shells  seem  to  respond  similarly  to  the  instantaneous  change  in  the  charge  on  the 
solute  and  this  response  follows  a  gaussian  behavior.  However,  the  first  shell  molecules  begin 
to  react  to  the  field  exerted  by  the  neighboring  dipoles  earlier  thAn  the  molecules  in  other 
shells.  Since  this  shell  contributes  more  than  50%  to  the  overall  response  fimction,  one 
expects  to  observe  the  characteristics  of  the  Ist  shell  relaxation  in  the  overall  response 
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function.  The  initial  relaxation  behavior  of  the  solvent  molecules  in  the  subsequent  shells 
does  not  show  any  appreciable  difference  from  each  other,  except  for  a  small  time  lag  in  the 
propagation  of  the  response.  These  shells  are  seen  to  be  completely  relaxed  within  1  ps  time 
period,  but  the  first  shell  remains  unrelaxed  for  another  picosecond.  The  emerging  picture 
is  quite  compatible  with  the  Qnsager  'snowball*  coiyecture. 

We  represent  a  aimilar  set  of  curves  obtained  fi^m  the  tnyectories  from  NESn  in 
Figure  6b.  Except  for  the  small  differences  observed  in  the  relaxation  of  the  first  shell,  the 
curves  are  very  similar  to  the  ones  in  Figure  6a.  Consequently,  the  relaxation  picture 
presented  above  is  also  valid  for  the  case  when  the  translational  motion  of  the  solvent 
participates  in  the  dynamics,  but  the  rotational  motion  is  dominant  (small  p’). 

In  Figure  6c,  we  plot  S(t)  of  the  solvent  shells  vs  time  for  the  case  in  which  solvent  is 
modified  to  have  T >2.5.  One  can  notice  the  difference  between  shells  even  in  the  behavior 
of  the  initial  fast  component  of  the  relaxation.  Thouj^  the  decay  is  gaussian,  the  curves  have 
a  wider  spread  compared  to  the  previous  cases.  All  the  curves  reach  their  equilibrium  values 
at  the  same  time  (around  7  ps),  thus  displaying  a  breakdown  of  the  (Dnsager  conjecture.  This 
agrees  with  the  predictions  of  Bagchi  and  Chandra  that  pointed  out  that  Oosager’s  conjecture 
becomes  incorrect  when  the  rotational  mode  is  not  a  dominant  mode  m  the  relaxation 
mechanism  [24]. 

For  a  confirmation  of  our  conclusions  let  us  now  turn  to  Figures  7.  In  these  we  display 
the  time  dependence  of  the  average  cosine  of  the  angle  (e)  between  the  solvent  dipole  and 
the  line  joining  the  centers  of  solute  and  sohwnt  We  can  easily  obtain  the  average  value  of 
cos(e)  for  a  particular  shell  at  any  time  t  Note,  that  for  a  neutral  solute  at  equilibrium  with 
the  solvent,  <co8(e)>  is  zero,  since  the  solvent  molecules  are  randomly  oriented  aroxmd  the 
solute.  However,  when  a  charge  is  placed  on  the  solute,  <cos(e)>  may  change  from  shell  to 
shell  and  the  value  of  <cos(e)>  for  a  particular  shell  depends  on  the  distance  from  the  ion. 

Figure  7a  represents  how  the  value  of  <cos(e)>  varies  with  time  for  different  solvent 
shells  when  the  translational  motion  is  absent.  At  t-0,  almost  all  the  shells  have  their 
<cos(e)>  values  dose  to  zero  showing  that  the  initial  configurations  of  the  non-equilibrium 
trajectories  are  selected  firom  a  rather  good  random  distribution,  ^th  the  time  proceeding, 
<co6(e)>  evolves  towards  values  obtained  fium  the  equilibrium  simulations. 

As  we  see  from  Figure  7a  the  molecules  in  the  first  shell  reorient  fast  towards  the  ion 
and  even  "overshoot”  the  equilibrium  angle.  Subsequently,  due  to  interactions  with  the 
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neighbors  and  the  ion  the  average  angle  oscillates  on  its  way  towards  equilibrium.  For  the 
molecules  in  other  shells  the  picture  is  similar,  but  due  to  the  smaller  perturbation  from  the 
ion  the  deviation  in  the  angle  is  smaller,  as  a  result  the  relaxation  proceeds  faster  to 
the  equilibrium.  Not  surprisingly,  the  values  of  the  tima  taken  to  complete  relaxation 
obtained  from  Figure  7a  agree  quite  well  with  those  evaluated  from  S(t)  plots. 

In  Figure  7b,  we  present  a  similar  set  of  plots  for  the  case  when  the  solvent 
transliitions  are  included  in  molecular  dynamics  (NESII).  Inclusion  of  the  translation  is  visible 
only  in  the  features  of  the  first  shell  and  the  final  (equilibrium)  values  of  the  first  three 
shells.  The  last  shell  displays  quite  the  same  results  for  the  two  cases.  The  change  in  the 
equilibrium  values  of  the  angle  is  due  to  the  electrostriction,  which  brings  shells  closer  to  the 
ion  and  the  distance  from  the  ion  determines  the  ani^e.  Observation  of  the  similar  behavior 
of  the  angular  relaxation  confirms  the  rotationally  dominant  nature  of  the  solvent  relaxation 
for  small  value  of  the  parameter  p*.  Somewhat  different  features  can  be  observed  in  the  plots 
of  <cos(e)>  vs  time  (Figure  7c)  for  the  case  in  which  T  a2.5.  For  this  case  all  the  shells  reach 
their  equilibrium  values  around  7p8  in  good  agreement  with  the  results  of  S(t)  plots.  Again 
it  proves  the  breakdown  of  the  Onsager  conjecture.  As  we  see  from  the  figure  the  first 
minimiim  of  <co6(0)>  does  not  deviate  much  firom  its  equilibrium  value  in  this  case.  This 
is  due  to  the  competition  between  translational  and  rotational  motions  which  contribute  in 
this  case  to  the  relaxation. 

Conchiaione 

We  studied  the  dynamics  of  ion  solvation  in  a  simple  Stockmayer  fluid  using  molecular 
dynamics  computer  simulation  technique.  From  our  simulations  we  have  concluded  that  the 
relaxation  proceeds  in  two  time  regimes:  initial  short  time  decay  regime,  that  can  be  called 
inertial  regime,  and  subsequent  long  time  decay  regime,  that  can  be  called  diffusional.  The 
short  tiwia  decay,  as  we  have  observed,  is  mostly  determined  by  the  relaxation  of  the  solvent 
molecules  belonging  to  the  first  solvation  shell  It  can  be  described  by  a  correlation  function, 
that  has  a  gaussian  form,  which  is  in  some  way  a  mathematical  reflection  of  the  inertial 
nature  of  the  motion  of  the  solvent  molecules.  We  have  also  observed  that  most  of  the  total 
energy  has  relaxed  during  this  short  time  deciQr.  The  inertial  character  of  the  initial  decay 
observed  in  our  qrstem  may  be  somewhat  dominant  because  of  the  nature  of  the  model 
desribing  the  solvent  (small  spherical  molecules  with  point  dipoles).  In  real  life  the  molecules 


are  not  spheres  and  the  liquids  are  more  hkely  to  be  in  the  overdamped  limit.  But  it  is  in 
some  way  surprising  that  all  the  simulations  up  to  date  show  the  dominanant  character  of  the 
inertial  motion  even  for  such  solvents  as  water  [11]. 

Our  simulations  also  show  that  the  loi^;  time  decay  of  the  relaxation  function  S(t)  has 
the  same  character  as  the  decay  of  the  correlation  function  of  the  single  dipole.  It  can, 
therefore,  be  characterized  by  an  exponential  relaxation  with  the  relaxation  time  r„  which 
is  known  to  be  of  the  same  order  as  [20].  This  is  why  we  think  of  this  time  regime  as  a 
difEUsional  regime.  Incidentally,  contrary  to  the  statements  in  the  literature,  we  observe  that 
the  difiEhsional  diaracter  of  the  long  time  decay  of  the  solvation  function  S(t)  is  not  related 
to  the  diffusional  rearrangement  of  the  first  solvation  shell  due  to  the  translational  difiusion 
of  the  solvent  molecules. 

We  have  also  investigated  the  validity  of  the  Onsager  conjecture  related  to  the 
"snowball”  efifect  We  have  observed  that  when  the  rotational  motion  is  the  only  mode  of  the 
relaxation  or  when  it  is  the  dominant  mode  of  the  relaxation  the  Onsager  coijecture  holds. 
When  the  role  of  translational  mode  is  increased  we  have  observed  that  the  Onsager 
conjecture  is  no  bnger  valid,  in  total  agreement  with  the  predictions  of  Bagchi  and  Chandra 
[10,241. 

We  observed  that  althou^  our  model  fluid  can  be  described  quite  well  by  a  simple 
Debye-type  continuum  model,  the  relaxation  process  can  not  be  described  by  a  continuum 
theory.  Even  more  sophisticated  theories  that  we  used  to  compare  with  our  results  fi*om  the 
simulatioiis,  (theories  based  on  dynamical  MSA  or  on  SVE)  are  not  performing  that  well 
either.  Very  recently  Bagchi  and  Chandra  presmited  a  non-Markovian  theory  which  included 
the  viscoelastic  behavior  of  the  solvent  responsible  for  the  initial  inertial  response  [26]. 
Numerical  studies  revealed  that  the  relaxation  function  displays  short-time  oscillations, 
followed  by  a  slow  long-time  decay  [26].  (Comparison  of  this  new  theory  with  the  computer 
simulation  performed  on  the  same  model  will  be  very  interesting. 

Note:  At  the  latest  stage  of  the  preparation  of  our  manuscript  we  learned  that  simulations 
on  similar  systems  were  performed  by  Prof.  Nitzan  and  his  collaborators  [27]. 
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Appoidix 

a.  The  functional  form  of  S(t)  in  the  framework  of  the  dynamical  MSA. 

The  dynamical  MSA  theory  was  proposed  by  Wol3mes  [6]  and  solved  exactly  by  Rips, 
{□after  and  Jortner  (RKJ)  [7]  and  by  Nichols  and  Calef  [8].  Here  we  present  the  final 
expression  for  the  relaxation  function  S(t)  we  use  in  order  to  get  the  corresponding  MSA 

curves.  For  the  detailes  on  the  development  of  the  MSA  the  reader  is  referred  to  the  original 
references  [6-8]. 

According  to  RKJ  we  can  write  an  analytical  e:q>re8sion  for  the  Laplace  transform  of 

A 

the  relaxation  function  S(p) 

S(p)  -  [&>)-xl[0)]/p[x(«o)-x(0)]  (Al) 


where  the  complex  admittance,  x(p)>  solvent  within  the  MSA  has  the  form: 


X<P)  -  Xmsa(p)  -  a-l/^))/2r,[UA(p)] 


(A2) 


The  dynamic  correction  factor  A(p)  is  given  ly  the  expression: 


A(p)  -  (3r,/r,){[^)]»/>+l/l!p)]-'/>-2}-' 


(A3) 


where 


(A4) 


and 


gi(p)  -  l+64[$p)]^^ 


(A5) 


rj  and  r,  are  the  solute  and  solvent  radii  respectively.  In  equations  (A2),  A(5)  the  dielectric 
susceptibility  of  the  medium  has  the  Debye  form: 

«(p)  -  ^.♦(q)'0/(l+pT4)  (A6) 
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Substitution  of  the  equations  (A2)-(A6)  into  equation  (Al)  results  in  the  functional  form  for 
the  relaxation  function  S(p).  To  get  the  desired  function  S(t)  the  inverse  Laplace  transform 
was  performed  with  the  help  of  numerical  techniques  [28]. 

b.  The  functional  form  of  S<t)  in  the  framework  of  SVE. 

The  Smoluchowski-Vlasov  Equation  (SVE)  was  originally  proposed  by  Wolynes  and 
Calef  [9]  to  describe  a  structural  relaxation  of  the  polarization  through  diffii«inn  in  the  mean 
field.  Subsequently  Nicols  and  Calef  [8]  solved  it  for  the  case  when  only  rotational  diffiiainn 
of  the  molecules  was  considered.  Bagchi  and  Chandra  [10b]  extended  the  solution  to  include 
the  translational  difiusion  of  the  molecules. 

The  fiinctional  form  of  S(t)  in  the  fimnework  of  SVE  is  given  by  the  equation  (19)  from 
ref.  10b,  which  is: 


m  m 

S(f)-2lirdjfeexp[-f/T^(*)]( 


(A7) 


where  ri^(k)  is  given  by  the  following  expression: 

l/ri(*)-2Da{l+p'(ife(7)»-l/3p[l+p'(;fea)»](C^+2CB)}  (A8) 


p  is  the  density  of  the  solvent,  is  the  solvent  rotational  difiusion  coefficient,  a  is  the 
solvent  diameter  and  p’  is  the  dimensionless  solvent  parameter  defined  by  equation  (7).  In 
equation  (A8)  C^  and  Cq  are  components  of  the  direct  correlation  function  C(k).  In  linear 
theories  C(k)  can  be  separated  into  two  parts,  Le. 


C(*)  -  CJ*CbD 


(A9) 
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where 


D  -  Zkk-I 


(AlO) 


If  MSA  is  used  to  find  C(k),  the  analytical  expressions  for  and  Cq  exist,  and  for  the  linear 
combination  that  appears  in  (A9)  we  get: 

C^*2Co  -  eKCp,(k,2Kp)  (All) 


where  Cpy(k;2Kp)  is  the  Percus-Yevick  direct  correlation  fimction  for  hard  spheres  at  density 
2Kp.  K  is  a  parameter  obtained  from  the  relationships: 

K  -  {/n  ^^12) 


and 


rj  -  irpo*/6 


(A13) 


where  (  is  the  solution  of  the  equation 


(l-25)» 


(A14) 


The  solution  of  this  equation  and  the  form  of  Ci>y0c;2Kp)  are  given  in  references  8  and  29. 
Note  that  in  our  numerical  work  we  used  r,«  r,  «  a/2. 
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Figure  Captions 

FIG.  1.  The  solute-solvent  radial  distribution  functions  obtained  from  ESI  and  ESn.  The 
solid  line  corresponds  to  the  ion-solvent  rdf^  the  dotted  line  to  the  neutral  solute- 
solvent  rdf. 

FIG.  2.  Orientational  distributions  for  different  solvation  shells  from  ESI  (upper  panel)  and 
ESn  Gower  panel);  Ist  shell  (-  2nd  shell  ( — ),  3rd  shell  (— )  and  bulk-like  (— ). 

FIG.  3.  Response  functions  S(t)  obtained  from  non-equilibrium  simulation  I  (NESI,  see  text) 
(— ),  dynamical  MSA  (-  •  SVER  ( — )  and  the  continuum  model  (—). 

FIG.  4.  a)  Response  functions  S(t)  obtained  from  NESI  ( — ),  NESn  (— )  and  SVERT  ( . ). 

b)  Response  function  obtained  from  NESITI. 

FIG.  5.  a)  Response  function  S(t)  [solid  line]  and  the  single  particle  dipole  autocorrelation 
function  *(t).  The  moment  of  inertia  in  the  simulations  is  T «0.025. 
b)  Same  as  above,  but  T «2.6. 

FIG.  6.  Response  functions  ( — )  and  their  solvent  shell  components  [1st  shell  (-•  2nd  shell 
( . ),  3rd  shell  (— )  and  bulk  (— )]  from  a)  NESI,  b)  NESII  and  c)  NESm. 

FIG.  7.  Variation  of  <cose(t)>  with  time  for  different  solvation  shells  [1st  shell  (— ), 

2nd  shell  (• — ),  3rd  shell  (— )  and  bulk  -)]  from  a)  NESI,  b)  NESII  and  c)  NESm. 


TABLE  1.  A  summary  of  mformation  on  solvation  shells.  The  values  in  the  table  are  obtained 
from  ESI  and  ESn  simnlatinns, 


a)  From  an  equilibrium  simulation  with  a  neutral  solute 


sheU 


radius  in 
reduced  units 


2 


3 


bulk-like 


number  of 
molecules 


1L8 


4L7 


80.6 


376.9 


<co8(e)> 

<density>  in 

reduced  units 

-0.0061 

0.798 

-0.0036 

0.800 

-0.0004 

0.825 

0.0003 

0.822 

b)  From  an  equilibrium  simulation  with  a  charged  solute 


solvation  shell  radius  in 

reduced  imits 


number  of 

nmlecules 

<cos(e)> 

9.5 

0.889 

33.3 

0.336 

76.4 

0.170 

39L8 

0.054 

<  density  >  in 
reduced  units 


1.226 


0.824 


0.817 


.816 
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